function [ SinTheta , WeightTheta ] = polar(NumPolar)
% function [ sin_theta_vect , WeightTheta ] = polar(NumPolar)
%
% This function returns the polar angle cosines and weights
% using the Tabuchi-Yamamoto quadrature set.
%
% References: 1) K. Smith, MIT NSE guest lecture notes
%
% J. Roberts,  02/27/2011
%
% Inputs:
%       NumPolar        - number of polar angles
%
% Outputs:
%       SinTheta        - vector of polar cosines
%       WeightTheta     - associated weights
%
SinTheta    = zeros(NumPolar,1);
WeightTheta = zeros(NumPolar,1);

% Assigning quadrature set
if NumPolar == 1
    SinTheta(1)         = 0.798184;
    WeightTheta(1)      = 1.0;
elseif NumPolar == 2
    SinTheta(1)         = 0.363900;
    SinTheta(2)         = 0.899900;
    WeightTheta(1)      = 0.212854;
    WeightTheta(2)      = 0.787146;
elseif NumPolar == 3
    SinTheta(1)         = 0.166648;
    SinTheta(2)         = 0.537707;
    SinTheta(3)         = 0.932954;
    WeightTheta(1)      = 0.046233;
    WeightTheta(2)      = 0.283619;
    WeightTheta(3)      = 0.670148;
elseif NumPolar == 8 % this uses gauss
    SinTheta = [	0.145209468829167
                    0.328295637781257
                    0.500682156115631
                    0.655258864753334
                    0.786275363090364
                    0.888943547917690
                    0.959530843792170
                    0.995476078554554 ];
    % 
    WeightTheta(1)      = 0.027152459411742137;
    WeightTheta(2)      = 0.06225352393877526;
    WeightTheta(3)      = 0.09515851168244673;
    WeightTheta(4)      = 0.12462897125537288;
    WeightTheta(5)      = 0.1495959888165501;
    WeightTheta(6)      = 0.16915651939499357;
    WeightTheta(7)      = 0.1826034150449233;
    WeightTheta(8)      = 0.1894506104550685;
elseif NumPolar == 16
   % i'm using the 1-d Sn contants, but solving using sqrt(1-mu^2) for the
   % sine of the polar angle
   SinTheta = [    0.073924216932330
        0.169026472215009
        0.263123146442321
        0.354895236900550
        0.443405441696658
        0.527801721048956
        0.607285351333596
        0.681108908323023
        0.748580189486052
        0.809067480925366
        0.862004920127778
        0.906897514663639
        0.943325622899806
        0.970948792813717
        0.989508894511019
        0.998832503193504 ];
    WeightTheta = [   0.007018600864948
        0.016274488699175
        0.025392015975798
        0.034273965082979
        0.042835820425217
        0.050998054624224
        0.058684061392098
        0.065822218819857
        0.072345791267224
        0.078193893915876
        0.083311921994799
        0.087652090660441
        0.091173876255960
        0.093844396569711
        0.095638717520179
        0.096540085931513 ];
  
else
    error(' ERROR: Polar quadrature for 1, 2, or 3 angles only! ');
end

end